#ifndef MATRIX_H
#define MATRIX_H
#include <iostream>
#include <sstream>
#include <stdlib.h>
#include <string>
#include <time.h>
#include <math.h>
#define pi 3.14
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-----Definição da classe para trabalhar com matrizes genéricas---//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
using namespace std;
template <class TipoClasse>
class matrix
{
public:
    int lin, col;
    TipoClasse** Mat;

public:
    matrix();
    matrix<TipoClasse>(const matrix<TipoClasse>& otherMatrix);
    ~matrix();
    // Acessando Atributos da classe
    int rows(){return this->lin;}
    int cols(){return this->col;}
    TipoClasse MatElem(int a, int b){return this->Mat[a][b];}
    //

    void init(int a, int b);
    void init(TipoClasse vetor[]);
    void init(TipoClasse *vetor, int LinCol);
    void init(TipoClasse **matriz, int a, int b);
    void init(string valor);
    void add(int a, int b, TipoClasse valor);
    void eye(int c);
    void zeros(int a, int b);
    void ones(int a, int b);
    void randU(int a, int b);//distribuição uniforme
    //sort
    TipoClasse max();
    TipoClasse min();

    //Operadores de Soma e Subtração
    matrix<TipoClasse> operator+ (matrix<TipoClasse> M1);//Soma matriz + matriz
    matrix<TipoClasse> operator+ (TipoClasse c);
    friend matrix<TipoClasse> operator+ (TipoClasse c, matrix <TipoClasse> M1){return M1+c;}
    matrix<TipoClasse> operator- (matrix<TipoClasse> M1);//Soma matriz + matriz
    matrix<TipoClasse> operator- (TipoClasse c);
    friend matrix<TipoClasse> operator- (TipoClasse c, matrix <TipoClasse> M1){return M1-c;}

    //Operadores de Multiplicação e Divisão
    matrix<TipoClasse> operator* (matrix<TipoClasse> M1);
    matrix<TipoClasse> operator* (TipoClasse c);
    friend matrix<TipoClasse> operator* (TipoClasse c,matrix <TipoClasse> M1){return M1*c;}
    matrix<TipoClasse> operator/ (matrix<TipoClasse> M1);//matriz vezes uma matriz invertida
    matrix<TipoClasse> operator/ (TipoClasse c);//divisão matriz escalar
    friend matrix<TipoClasse> operator/ (TipoClasse c,matrix <TipoClasse> M1){return c*(M1^-1);}

    //Potência de Matrizes
    matrix<TipoClasse> operator^ (int c);
    matrix<TipoClasse> operator^= (int c);
    //Operadores de Concatenação
    matrix<TipoClasse> operator&(matrix<TipoClasse> M1);//concatenação de linhas
    matrix<TipoClasse> operator|(matrix<TipoClasse> M1);//concatenação de colunas

    void operator= (TipoClasse *vetor);
    void operator= (matrix<TipoClasse> M1);
    void operator= (string valor);
    matrix<TipoClasse> operator~();//transposta de uma matriz
    // Álgebra Linear

    TipoClasse det();//determinante de uma matriz
    matrix<TipoClasse> transpose();//transposta de uma matriz
    matrix<TipoClasse> cofactor();//Matrizde cofatores: obs a adjunta é sua transposta
    matrix<TipoClasse> inv();//Inversa Matricial
    TipoClasse trace();//Traço da matriz
    matrix<TipoClasse> Gauss(matrix<TipoClasse> B);//Método da Eliminação de Gauss
    matrix<TipoClasse> LU();//Fatoração LU
    matrix<TipoClasse> eig();//Autovalores de uma matriz

    //Seleção de linhas e colunas de matrizes
    matrix<TipoClasse> selec(int a, string b);
    matrix<TipoClasse> selec(string a, int b);
    matrix<TipoClasse> selec(string a, string b);
    void selec(string a, string b, matrix<TipoClasse> M1, string c, string d);
    void selec(int a, string b, matrix<TipoClasse> M1, string c, string d);
    void selec(string a, int b, matrix<TipoClasse> M1, string c, string d);
    void selec(string a, string b, matrix<TipoClasse> M1, int c, string d);
    void selec(string a, string b, matrix<TipoClasse> M1, string c, int d);
    void selec(int a, int b, matrix<TipoClasse> M1, string c, string d);
    void selec(int a, string b, matrix<TipoClasse> M1, int c, string d);
    void selec(int a, string b, matrix<TipoClasse> M1, string c, int d);
    void selec(string a, int b, matrix<TipoClasse> M1, int c, string d);
    void selec(string a, int b, matrix<TipoClasse> M1, string c, int d);
    void selec(string a, string b, matrix<TipoClasse> M1, int c, int d);
    void selec(string a, int b, matrix<TipoClasse> M1, int c, int d);
    void selec(int a, string b, matrix<TipoClasse> M1, int c, int d);
    void selec(int a, int b, matrix<TipoClasse> M1, string c, int d);
    void selec(int a, int b, matrix<TipoClasse> M1, int c, string d);
    void selec(int a, int b, matrix<TipoClasse> M1, int c, int d);

    //-----------------------------------------------------------------//
    //-----------------------------------------------------------------//
    //------------Operações sob os elementos das matrizes--------------//
    //-----------------------------------------------------------------//
    //-----------------------------------------------------------------//

    friend matrix<TipoClasse> cos(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = cos(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> sin(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = sin(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> tan(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = tan(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> acos(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = acos(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> asin(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = asin(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> atan(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = atan(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> atan2(matrix<TipoClasse> M1, matrix<TipoClasse> M2)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = atan2(M1.Mat[i][j], M2.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> cosh(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = cosh(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> sinh(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = sinh(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> tanh(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = tanh(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> acosh(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = acosh(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> asinh(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = asinh(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> atanh(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = atanh(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> exp(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = exp(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> log(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = log(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> log10(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = log10(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> ceil(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = ceil(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> floor(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = floor(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> round(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col; j++)
                ret.Mat[i][j] = round(M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> diff(matrix<TipoClasse> M1)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col - 1; j++)
                ret.Mat[i][j] = (M1.Mat[i][j+1] - M1.Mat[i][j]);

        return ret;

    }

    friend matrix<TipoClasse> diff(matrix<TipoClasse> M1, TipoClasse h)
    {
        matrix<TipoClasse> ret = M1;
        for(  int i = 0; i < M1.lin; i++)
            for(  int j = 0; j < M1.col - 1; j++)
                ret.Mat[i][j] = (M1.Mat[i][j+1] - M1.Mat[i][j])/h;

        return ret;

    }

    void print();
};//Fim da classe;



//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-----------------Construto e Destrutor da classe-----------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//

template <class TipoClasse>
matrix<TipoClasse>::matrix()
{
    this->lin = 0;
    this->col = 0;
}

template <class TipoClasse>
matrix<TipoClasse>::matrix(const matrix<TipoClasse>& otherMatrix)
{
    this->init( otherMatrix.lin,otherMatrix.col);

    for(  int i = 0; i < this->lin; i++)
        for(  int j = 0; j < this->col; j++)
            this->Mat[i][j] = otherMatrix.Mat[i][j];

}

template <class TipoClasse>
matrix<TipoClasse>::~matrix()
{

    try
    {
        for(int i = 0; i < this->lin; i++)
        {
            if(this->Mat[i]!= NULL)
            {
                free(this->Mat[i]);
                this->Mat[i] = NULL;
            }
        }

        free(this->Mat);
        this->lin = 0;
        this->col = 0;
        this->Mat = NULL;
    }
    catch(int e)
    {
         cout << "Matriz ja foi deletada";
    }

}

//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-------------------Inicializações das matrizes-------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//


template <class TipoClasse>
void matrix<TipoClasse>::init( int a, int  b)
{
      int i;

      if(a!=0 && b!=0)
      {
    //allocate the 2d array

        this->Mat = (TipoClasse**)malloc((a) * sizeof(TipoClasse*));

        for(i = 0; i < a; i++)
          this->Mat[i] = (TipoClasse*)malloc((b) * sizeof(TipoClasse));
      }
      else
      {
          cout << "Erro de alocação da matriz\n";
      }

      this->lin = a;
      this->col = b;
}

template <class TipoClasse>
void matrix<TipoClasse>::init(TipoClasse **matriz, int a,   int b)
{
    int i;
    int j;
    this->init(a,b);

    for(i = 0; i < a; i++)
        for(j = 0; j < b; j++)
            this->Mat[i][j] = matriz[i][j];
}

template <class TipoClasse>
void matrix<TipoClasse>::init(string valor)
{
    int contRow = 1, contCol = 1, positionComma = 0, positionSemicomma = 0, stopWhile = valor.length();
    string temp;

    while(contCol + contRow <= stopWhile)
    {
        positionComma = valor.find(",");
        positionSemicomma = valor.find(";");

        if(((positionComma < positionSemicomma) && (positionComma != -1)) || ((positionComma > positionSemicomma) && (positionSemicomma == -1)) || ((positionComma == -1) && (positionSemicomma == -1) && valor.length() != 0))
        {
            if(positionComma == -1 && valor.length())
            {
                positionComma = valor.length();
            }
            temp = valor.substr(0,positionComma);
            TipoClasse Itemp = (TipoClasse) atof(temp.c_str());
            this->add(contRow, contCol, Itemp);
            valor.erase(0, positionComma + 1);
            contCol++;
        }
        else if((positionComma > positionSemicomma) || ((positionComma < positionSemicomma) && (positionComma == -1)))
        {
            temp = valor.substr(0,positionSemicomma);
            TipoClasse Itemp = (TipoClasse) atof(temp.c_str());
            this->add(contRow, contCol, Itemp);
            valor.erase(0, positionSemicomma + 1);
            contCol = 1;
            contRow++;
        }
        else
            break;

    }

}

template <class TipoClasse>
void matrix<TipoClasse>::init(TipoClasse vetor[])
{
    int a = sizeof(vetor) + 1;
    this->init(a,1);
    for(int i = 0; i < this->lin; i++)
    {
        TipoClasse *p1 = &vetor[i];
        this->Mat[i][0] = p1[0];
    }
}

template <class TipoClasse>
void matrix<TipoClasse>::operator= (TipoClasse *vetor)
{
    this->init(vetor);
}

template <class TipoClasse>
void matrix<TipoClasse>::init(TipoClasse *vetor, int LinCol)
{
    int a = sizeof(vetor) + 1;
    if(LinCol == 1)
    {
        this->init(a,1);
        for(int i = 0; i < this->lin; i++)
        {
            TipoClasse *p1 = &vetor[i];
            this->Mat[i][0] = p1[0];
        }
    }
    else if(LinCol == 2)
    {
        this->init(1,a);
        for(int i = 0; i < this->col; i++)
        {
            TipoClasse *p1 = &vetor[i];
            this->Mat[0][i] = p1[0];
        }
    }
    else
        cout << "Erro de Inicialização do vetor";
}

template <class TipoClasse>
void matrix<TipoClasse>::add(  int a,   int b, TipoClasse valor)
{
    matrix<TipoClasse> A;
    if(a >= this->lin && b >= this->col)
        A.init(a,b);
    else if (a <= this->lin && b >= this->col)
        A.init(this->lin,b);
    else if (a >= this->lin && b <= this->col)
        A.init(a,this->col);
    else
        A.init(this->lin,this->col);

    for(  int i = 0; i < A.lin; i++)
        for(  int j = 0; j < A.col; j++)
            A.Mat[i][j] = 0;

    for(  int i = 0; i < this->lin; i++)
        for(  int j = 0; j < this->col; j++)
            A.Mat[i][j] = this->Mat[i][j];

    A.Mat[a-1][b-1] = valor;

    this->init(A.lin,A.col);

    for(  int i = 0; i < this->lin; i++)
        for(  int j = 0; j < this->col; j++)
            this->Mat[i][j] = A.Mat[i][j];

}


template <class TipoClasse>
void matrix<TipoClasse>::eye(  int c)
{
  this->init(c, c);

  for(unsigned int i = 0; i < this->lin; i++)
    for(unsigned int j = 0; j < this->col; j++)
      if (i == j)
        this->Mat[i][j] = 1;
      else
         this->Mat[i][j] = 0;
}

template <class TipoClasse>
void matrix<TipoClasse>::zeros(  int a, int b)
{
  this->init(a, b);

  for(unsigned int i = 0; i < this->lin; i++)
    for(unsigned int j = 0; j < this->col; j++)
        this->Mat[i][j] = 0;
}

template <class TipoClasse>
void matrix<TipoClasse>::ones(  int a, int b)
{
    this->init(a, b);

    for(unsigned int i = 0; i < this->lin; i++)
        for(unsigned int j = 0; j < this->col; j++)
            this->Mat[i][j] = 1;
}

template <class TipoClasse>
void matrix<TipoClasse>::randU(int a, int b)
{
    double t;
    srand (time(NULL));
    this->init(a, b);

    for(unsigned int i = 0; i < this->lin; i++)
        for(unsigned int j = 0; j < this->col; j++)
        {
            t=rand();
            this->Mat[i][j] = (TipoClasse)(t)/RAND_MAX;
        }
}


template <class TipoClasse>
TipoClasse matrix<TipoClasse>::max()
{
  unsigned int i, j;
  TipoClasse  maximum = this->Mat[0][0];

  for(i = 0; i < this->lin; i++)
    for(j = 0; j < this->col; j++)
      if(maximum < this->Mat[i][j])
        maximum = this->Mat[i][j];

  return maximum;
}

template <class TipoClasse>
TipoClasse matrix<TipoClasse>::min()
{
  unsigned int i, j;
  double  minimun;
  minimun = this->Mat[0][0];

  for(i = 0; i < this->lin; i++)
    for(j = 0; j < this->col; j++)
      if(minimun > this->Mat[i][j])
        minimun = this->Mat[i][j];

  return minimun;
}

//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-----------------Operadores de Soma e Subtração------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>:: operator+(matrix<TipoClasse> M1)
{
    matrix<TipoClasse> ret;
    if(this->lin == M1.lin && this->col == M1.col)
    {
        ret.init(this->lin,this->col);

        for(  int i = 0; i < this->lin; i++)
            for(  int j = 0; j < this->col; j++)
                 ret.Mat[i][j]= this->Mat[i][j] + M1.Mat[i][j];
    }
    else
    {
         cout << "Matrizes de dimensões diferentes";
    }

    return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator+(TipoClasse c)
{
    matrix<TipoClasse> ret;
    ret.init(this->lin, this->col);

  for(int i = 0; i < this->lin; i++)
    for(unsigned int j = 0; j < this->col; j++)
      ret.Mat[i][j] = this->Mat[i][j] + c;

  return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>:: operator-(matrix<TipoClasse> M1)
{
    matrix<TipoClasse> ret;
    if(this->lin == M1.lin && this->col == M1.col)
    {
        ret.init(this->lin,this->col);

        for(  int i = 0; i < this->lin; i++)
            for(  int j = 0; j < this->col; j++)
                 ret.Mat[i][j]= this->Mat[i][j] - M1.Mat[i][j];
    }
    else
    {
         cout << "Erro no comando";
    }

    return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator-(TipoClasse c)
{
    matrix<TipoClasse> ret;
    ret.init(this->lin, this->col);

  for(unsigned int i = 0; i < this->lin; i++)
    for(unsigned int j = 0; j < this->col; j++)
      ret.Mat[i][j] = this->Mat[i][j] - c;

  return ret;
}

//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-------------Operadores de Multiplicação e Divisão---------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator*(matrix<TipoClasse> M1)
{
    matrix<TipoClasse> ret;
    if(this->col == M1.lin)
    {
        ret.init(this->lin,M1.col);

        for(  int i = 0; i < this->lin; i++)
          for(  int j = 0; j < M1.col; j++)
            for(  int inner = 0; inner < M1.lin; inner++)
                ret.Mat[i][j] = (TipoClasse) (ret.Mat[i][j] + (this->Mat[i][inner] * M1.Mat[inner][j]));

    }
    else
    {
         cout << "Erro no comando";
    }

    return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator*(TipoClasse c)
{
    matrix ret;
    ret.init(this->lin, this->col);

  for(unsigned int i = 0; i < this->lin; i++)
    for(unsigned int j = 0; j < this->col; j++)
      ret.Mat[i][j] = (TipoClasse) (this->Mat[i][j] * c);

  return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator/(matrix<TipoClasse> M1)
{
    matrix<TipoClasse> ret;
    ret.init(this->lin, this->col);

    for(int i = 0; i < this->lin; i++)
        for(int j = 0;j < this->col; j++)
            ret.Mat[i][j] = this->Mat[i][j];

    ret = ret*M1.inv();

    return ret;

}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator/(TipoClasse c)
{
    matrix ret;
    ret.init(this->lin, this->col);

  for(unsigned int i = 0; i < this->lin; i++)
    for(unsigned int j = 0; j < this->col; j++)
      ret.Mat[i][j] = (TipoClasse) (this->Mat[i][j] / c);

  return ret;
}

//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//---------------------- Potência de Matrizes ---------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator^(int c)
{
  matrix<TipoClasse> ret, Temp;
  Temp.init(this->lin, this->col);
  ret.eye(this->lin);

  for(int i = 0; i < this->lin; i++)
    for(int j = 0; j < this->col; j++)
      Temp.Mat[i][j] = this->Mat[i][j];

  if(this->lin != this->col)
    return ret;

  unsigned int n = this->lin;

  if(c > 0)
    for(unsigned int i=0; i<c; i++)
      ret = ret*Temp;
  else if(c == 0) {
    ret.eye(n);
    return ret;
  }
  else {
    c = -c;
    for(unsigned int i=0; i<c; i++)
      ret = ret*Temp;

    ret = ret.inv();
  }

    return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator^=(int c)
{
  matrix<TipoClasse> ret;
  ret.init(this->lin, this->col);

  for(int i = 0; i < this->lin; i++)
    for(int j = 0; j < this->col; j++)
      ret.Mat[i][j] = pow(this->Mat[i][j],c);

    return ret;
}
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//------------------- Operadores de Concatenação ------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//

//ConcatenaÃ§Ã£o em linhas
template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator&(matrix<TipoClasse> M1)
{
    matrix<TipoClasse> ret;
    int temp = this->col + M1.col;

    if(this->lin == M1.lin)
    {
        ret.init(this->lin, temp);
        for(unsigned int i = 0; i < this->lin; i++)
        for(unsigned int j = 0; j < temp; j++)
        {
            if(j < this->col)
                ret.Mat[i][j] = this->Mat[i][j];
            else
                ret.Mat[i][j] = M1.Mat[i][j-this->col];
        }
    }
    else if(this->lin == 0 && this->col == 0)
        ret = M1;
    else
    {
    //  Essa parte precisa ser melhor elaborada
    }

    return ret;
}

//ConcatenaÃ§Ã£o em colunas
template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator|(matrix<TipoClasse> M1)
{
  matrix<TipoClasse> ret;
  unsigned int temp = this->lin + M1.lin;
  ret.init(temp, this->col);

  if(this->col == M1.col)
    for(unsigned int i = 0; i < temp; i++)
      for(unsigned int j = 0; j < ret.col; j++) {
        if(i < this->lin)
          ret.Mat[i][j] = this->Mat[i][j];
        else
          ret.Mat[i][j] = M1.Mat[i-this->lin][j];
      }
  else if(this->lin == 0 && this->col == 0)
    ret = M1;
  else {
//  Essa parte precisa ser melhor elaborada
  }

  return ret;
}


//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//------------------------- Álgebra Linear ------------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//



template <class TipoClasse>
TipoClasse matrix<TipoClasse>::det()
{
    unsigned int i, j, j1, j2;
    double det = 0;
    matrix<TipoClasse> Temp;
    unsigned n = this->lin;

    if(this->lin != this->col)
        return (TipoClasse) 0.0;

    if (n < 1)
    { /* Error */

    }
    else if (n == 1)
    { /* Shouldn't get used */
        det = this->Mat[0][0];
    }
    else if (n == 2)
    {
        det = this->Mat[0][0] * this->Mat[1][1] - this->Mat[1][0] * this->Mat[0][1];
    }
    else
    {
        for(j1=0; j1<n; j1++)
        {
          Temp.zeros(n-1,n-1);
          for(i=1; i<n; i++)
          {
            j2 = 0;
            for(j=0; j<n; j++)
            {
              if(j == j1)
                continue;

              Temp.Mat[i-1][j2] = this->Mat[i][j];
              j2++;
            }
          }
          det += (TipoClasse) pow(-1.0, j1+2.0) * this->Mat[0][j1] * Temp.det();
    }
    }

    return(det);

}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::transpose()
{
  matrix<TipoClasse> trans;
  trans.zeros(this->col, this->lin);

  for(int i = 0; i < this->lin; i++)
    for(int j = 0; j < this->col; j++)
      trans.Mat[j][i] = this->Mat[i][j];

  return trans;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::operator~()
{
    matrix<TipoClasse> trans;
    trans.zeros(this->col, this->lin);

    for(int i = 0; i < this->lin; i++)
      for(int j = 0; j < this->col; j++)
        trans.Mat[j][i] = this->Mat[i][j];

    return trans;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::cofactor()
{
  unsigned int i, j, ii, jj, i1, j1;
  TipoClasse det;
  matrix<TipoClasse> Temp, ret;
  unsigned n = this->lin ;
  Temp.zeros(n-1, n-1);
  ret.zeros(n, n);

  if(this->lin != this->col)
    return ret;
  if (n < 1) { /* Error */

  }

  for(j=0; j<n; j++)
    for(i=0; i<n; i++) {

      /* Form the adjoint a_ij */
      i1 = 0;
      for(ii=0; ii<n; ii++) {
        if(ii == i)
          continue;

        j1 = 0;
        for(jj=0; jj<n; jj++) {
          if(jj == j)
            continue;

          Temp.Mat[i1][j1] = this->Mat[ii][jj];
          j1++;
        }
        i1++;
      }

      /* Fill in the elements of the cofactor */
      ret.Mat[i][j] = (TipoClasse) pow(-1.0, i+j+2.0) * Temp.det();
    }

  return ret;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::inv()
{
  matrix<TipoClasse> ret;
  ret = this->cofactor();
  ret = ret.transpose() / this->det();

  return ret;
}

template <class TipoClasse>
TipoClasse matrix<TipoClasse>::trace()
{
    TipoClasse tr = 0;

    if(this->lin != this->col)
      return tr;

    for(unsigned int i=0; i < this->lin; i++)
      tr += this->Mat[i][i];

    return tr;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::Gauss(matrix<TipoClasse> B)
{
  unsigned int i, j, k;
  matrix<TipoClasse> Temp, ret, X;
  Temp.zeros(this->lin, this->col);
  ret.zeros(this->lin,this->col+1);

  for(i = 0; i < this->lin; i++)
    for(j = 0; j < this->col; j++)
      Temp.Mat[i][j] = this->Mat[i][j];

  ret = Temp&B;
//  ret.print();
  for(j = 0; j < ret.lin; j++)
  {
      if(ret.Mat[j][j] == 0 && j < ret.lin - 1)
      {
          matrix<TipoClasse> Ttemp = ret.selec(j+1,":");
//          Ttemp.print();
          ret.selec(j+1,":",ret,j+2,":");
//          ret.print();
          ret.selec(j+2,":",Ttemp,":",":");// a função seleciona a linha j+2 de ret e coloca no lugar dela todas as linhas e colunas de Ttemp
//          cout << endl;
//          ret.print();
      }
      if(ret.Mat[j][j] == 0 && j == ret.lin - 1)
      {
          matrix<TipoClasse> Ttemp = ret.selec(j+1,":");
          ret.selec(j+1,":",ret,j,":");
          Ttemp.selec(j,":",ret,":",":");
      }

  }
//  ret.print();

//  cout<< ret.lin<<' '<<ret.col<<'\n';
  for(j = 0; j < ret.col-2; j++)
      for(i = j; i < ret.lin-1; i++) {
          double m = ret.Mat[i+1][j]/ret.Mat[j][j];
          for(k = j; k < ret.col; k++)
            ret.Mat[i+1][k] = ret.Mat[i+1][k] - m*ret.Mat[j][k];
      }

  X.zeros(B.lin,B.col);
  double R;
  for(i = 0; i < ret.lin; i++) {
    int k = ret.lin - i - 1;
    R=0;
    for(j = k+1; j < ret.col-1; j++)
      R = R + ret.Mat[k][j]*X.Mat[j][0];

    X.Mat[k][0]=(ret.Mat[k][ret.col-1]-R)/ret.Mat[k][k];
  }
  return X;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::LU()
{
  unsigned int i, j, k;
  matrix<TipoClasse> U, L, LU;
  U.zeros(this->lin, this->col);
  L.eye(this->lin);
  LU.zeros(this->lin,this->col+this->col);

  for(i = 0; i < this->lin; i++)
    for(j = 0; j < this->col; j++)
      U.Mat[i][j] = this->Mat[i][j];

  for(j = 0; j < U.lin; j++)
  {
      if(U.Mat[j][j] == 0 && j < U.lin - 1)
      {
          matrix<TipoClasse> Ttemp = U.selec(j+1,":");
//          Ttemp.print();
          U.selec(j+1,":",U,j+2,":");
//          ret.print();
          U.selec(j+2,":",Ttemp,":",":");// a função seleciona a linha j+2 de U e coloca no lugar dela todas as linhas e colunas de Ttemp
//          cout << endl;
//          ret.print();
      }
      if(U.Mat[j][j] == 0 && j == U.lin - 1)
      {
          matrix<TipoClasse> Ttemp = U.selec(j+1,":");
          U.selec(j+1,":",U,j,":");
          Ttemp.selec(j,":",U,":",":");
      }

  }

//   cout<< ret.rows<<' '<<ret.col<<'\n';
  for(j= 0; j < U.col-1; j++)
      for(i = j; i < U.lin-1; i++) {
          double m = U.Mat[i+1][j]/U.Mat[j][j];
          for(k = j; k < U.col; k++){
//             cout<<i+1<<' '<<j<<' '<<k<<'\n';
            U.Mat[i+1][k] = U.Mat[i+1][k] - m*U.Mat[j][k];
          }
          L.Mat[i+1][j] = m;
      }
  LU = L,U;
  return LU;
}

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::eig()
{
  matrix<TipoClasse> Aav, TempLU, L, U, eigv;
  U.zeros(this->lin, this->col);
  L.zeros(this->lin, this->col);
  Aav.zeros(this->lin, this->col);
  TempLU.zeros(this->lin,this->col+this->col);
  eigv.zeros(this->lin, 1);

  double ia,ip, it;
  unsigned int i, j;

  for(i = 0; i < this->lin; i++)
    for(j = 0; j < this->col; j++)
      Aav.Mat[i][j] = this->Mat[i][j];

  i = 0;
  while(true){
    ia = Aav.max();
    TempLU = Aav.LU();
    for(i = 0; i < TempLU.lin; i++)
      for(j = 0; j < TempLU.col; j++) {
        if(j < this->col)
          L.Mat[i][j] = TempLU.Mat[i][j];
        else
          U.Mat[i][j - this->col] = TempLU.Mat[i][j];
      }
    Aav = U*L;
    ip = Aav.max();
    it = (ia-ip);
    if(it < 0)
      it = -it;
    if ( it < 0.0001 || i >= 1000)
        break;
    i++;
  }
  for(i = 0; i < this->lin; i++)
    eigv.Mat[i][0] = Aav.Mat[i][i];

  return eigv;
}

//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-------------------Seleção de Linhas e Colunas-------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//

//Seleciona de uma linha um numero específicos de colunas

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::selec(int a, string b)
{
    matrix<TipoClasse> ret;
    string temp;
    int position = b.find(":");
    int iniLin = a - 1, finalLin = a - 1, iniCol = 0, finalCol = this->col-1;

//correção ainda em teste{
    if(iniLin < 0)
            iniLin = 0;

    if(finalLin < 0)
            finalLin = 0;

    if(finalCol < 0)
            finalCol = 0;
//correção ainda em teste }

    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (colunas)
        temp = b.substr(0,position);
        iniCol = atoi(temp.c_str()) - 1;
        temp = b.substr(position + 1, b.length() - position + 1);
        finalCol = atoi(temp.c_str()) - 1;
    }


    //Selecionando a submatriz resultante

    ret.init(finalLin - iniLin + 1, finalCol - iniCol + 1);

    for(int i = iniLin; i <= finalLin; i++)
        for(int j = iniCol; j <= finalCol; j++)
            ret.Mat[i - iniLin][j - iniCol] = this->Mat[i][j];

    return ret;
}

//Seleciona de uma coluna um numero específicos de linhas

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::selec(string a, int b)
{
    matrix<TipoClasse> ret;
    string temp;
    int position = a.find(":");
    int iniLin = 0, finalLin = this->lin-1, iniCol = b - 1, finalCol = b - 1;
    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (linhas)
        temp = a.substr(0,position);
        iniLin = atoi(temp.c_str()) - 1;
        temp = a.substr(position + 1, a.length() - position + 1);
        finalLin = atoi(temp.c_str()) - 1;
    }

    //Selecionando a submatriz resultante

    ret.init(finalLin - iniLin + 1, finalCol - iniCol + 1);

    for(int i = iniLin; i <= finalLin; i++)
        for(int j = iniCol; j <= finalCol; j++)
            ret.Mat[i - iniLin][j - iniCol] = this->Mat[i][j];

    return ret;
}

//Seleciona um numero específico de linhas e colunas de uma matriz

template <class TipoClasse>
matrix<TipoClasse> matrix<TipoClasse>::selec(string a, string b)
{
    matrix<TipoClasse> ret;
    string temp;
    int position = a.find(":");
    int iniLin = 0, finalLin = this->lin-1, iniCol = 0, finalCol = this->col-1;
    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (linhas)
        temp = a.substr(0,position);
        iniLin = atoi(temp.c_str()) - 1;
        temp = a.substr(position + 1, a.length() - position + 1);
        finalLin = atoi(temp.c_str()) - 1;
    }

    position = b.find(":");
    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (colunas)
        temp = b.substr(0,position);
        iniCol = atoi(temp.c_str()) - 1;
        temp = b.substr(position + 1, b.length() - position + 1);
        finalCol = atoi(temp.c_str()) - 1;
    }


    //Selecionando a submatriz resultante

    ret.init(finalLin - iniLin + 1, finalCol - iniCol + 1);

    for(int i = iniLin; i <= finalLin; i++)
        for(int j = iniCol; j <= finalCol; j++)
            ret.Mat[i - iniLin][j - iniCol] = this->Mat[i][j];

    return ret;
}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, string b, matrix<TipoClasse> M1, string c, string d)
{
    matrix<TipoClasse> ret = M1;
    string temp;
    int position = a.find(":");
    int iniLinA = 0, finalLinA = this->lin-1, iniColA = 0, finalColA = this->col-1;
    int iniLinB = 0, finalLinB = M1.lin-1, iniColB = 0, finalColB = M1.col-1;

    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (linhas)
        temp = a.substr(0,position);
        iniLinA = atoi(temp.c_str()) - 1;
        temp = a.substr(position + 1, a.length() - position + 1);
        finalLinA = atoi(temp.c_str()) - 1;
    }

    position = b.find(":");
    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (colunas)
        temp = b.substr(0,position);
        iniColA = atoi(temp.c_str()) - 1;
        temp = b.substr(position + 1, b.length() - position + 1);
        finalColA = atoi(temp.c_str()) - 1;
    }

    position = c.find(":");

    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (linhas)
        temp = c.substr(0,position);
        iniLinB = atoi(temp.c_str()) - 1;
        temp = c.substr(position + 1, c.length() - position + 1);
        finalLinB = atoi(temp.c_str()) - 1;
    }

    position = d.find(":");
    if(position != 0)
    {
    //Encontrando os limites da submatriz resultante (colunas)
        temp = d.substr(0,position);
        iniColB = atoi(temp.c_str()) - 1;
        temp = d.substr(position + 1, d.length() - position + 1);
        finalColB = atoi(temp.c_str()) - 1;
    }


    for(int i = 0; i <= finalLinA - iniLinA ; i++)
        for(int j = 0; j <= finalColA - iniColA; j++)
            this->Mat[i + iniLinA][j + iniColA] = ret.Mat[i + iniLinB][j + iniColB];


}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, string b, matrix<TipoClasse> M1, string c, string d)
{
    stringstream temp;
    string str;

    temp << a;
    temp >> str;

    str = str+":"+str;
    this->selec(str, b, M1, c, d);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, int b, matrix<TipoClasse> M1, string c, string d)
{
    stringstream temp;
    string str;

    temp << b;
    temp >> str;

    str = str+":"+str;
    this->selec(a, str, M1, c, d);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, string b, matrix<TipoClasse> M1, int c, string d)
{
    stringstream temp;
    string str;

    temp << c;
    temp >> str;

    str = str+":"+str;
    this->selec(a, b, M1, str, d);

}


template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, string b, matrix<TipoClasse> M1, string c, int d)
{
    stringstream temp;
    string str;

    temp << d;
    temp >> str;

    str = str+":"+str;
    this->selec(a, b, M1, c, str);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, int b, matrix<TipoClasse> M1, string c, string d)
{
    stringstream temp1, temp2;
    string str1, str2;

    temp1 << a;
    temp1 >> str1;

    temp2 << b;
    temp2 >> str2;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    this->selec(str1, str2, M1, c, d);

}


template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, string b, matrix<TipoClasse> M1, int c, string d)
{
    stringstream temp1, temp2;
    string str1, str2;

    temp1 << a;
    temp1 >> str1;

    temp2 << c;
    temp2 >> str2;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    this->selec(str1, b, M1, str2, d);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, string b, matrix<TipoClasse> M1, string c, int d)
{
    stringstream temp1, temp2;
    string str1, str2;

    temp1 << a;
    temp1 >> str1;

    temp2 << d;
    temp2 >> str2;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    this->selec(str1, b, M1, c, str2);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, int b, matrix<TipoClasse> M1, int c, string d)
{
    stringstream temp1, temp2;
    string str1, str2;

    temp1 << b;
    temp1 >> str1;

    temp2 << c;
    temp2 >> str2;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    this->selec(a, str1, M1, str2, d);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, int b, matrix<TipoClasse> M1, string c, int d)
{
    stringstream temp1, temp2;
    string str1, str2;

    temp1 << b;
    temp1 >> str1;

    temp2 << d;
    temp2 >> str2;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    this->selec(a, str1, M1, c, str2);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, string b, matrix<TipoClasse> M1, int c, int d)
{
    stringstream temp1, temp2;
    string str1, str2;

    temp1 << c;
    temp1 >> str1;

    temp2 << d;
    temp2 >> str2;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    this->selec(a, b, M1, str1, str2);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(string a, int b, matrix<TipoClasse> M1, int c, int d)
{
    stringstream temp1, temp2, temp3;
    string str1, str2, str3;

    temp1 << b;
    temp1 >> str1;

    temp2 << c;
    temp2 >> str2;

    temp3 << d;
    temp3 >> str3;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    str3 = str3+":"+str3;
    this->selec(a, str1, M1, str2, str3);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, string b, matrix<TipoClasse> M1, int c, int d)
{
    stringstream temp1, temp2, temp3;
    string str1, str2, str3;

    temp1 << a;
    temp1 >> str1;

    temp2 << c;
    temp2 >> str2;

    temp3 << d;
    temp3 >> str3;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    str3 = str3+":"+str3;
    this->selec(str1, b, M1, str2, str3);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, int b, matrix<TipoClasse> M1, string c, int d)
{
    stringstream temp1, temp2, temp3;
    string str1, str2, str3;

    temp1 << a;
    temp1 >> str1;

    temp2 << b;
    temp2 >> str2;

    temp3 << d;
    temp3 >> str3;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    str3 = str3+":"+str3;
    this->selec(str1, str2, M1, c, str3);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, int b, matrix<TipoClasse> M1, int c, string d)
{
    stringstream temp1, temp2, temp3;
    string str1, str2, str3;

    temp1 << a;
    temp1 >> str1;

    temp2 << b;
    temp2 >> str2;

    temp3 << c;
    temp3 >> str3;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    str3 = str3+":"+str3;
    this->selec(str1, str2, M1, str3, d);

}

template <class TipoClasse>
void matrix<TipoClasse>::selec(int a, int b, matrix<TipoClasse> M1, int c, int d)
{
    stringstream temp1, temp2, temp3, temp4;
    string str1, str2, str3, str4;

    temp1 << a;
    temp1 >> str1;

    temp2 << b;
    temp2 >> str2;

    temp3 << c;
    temp3 >> str3;

    temp4 << d;
    temp4 >> str4;

    str1 = str1+":"+str1;
    str2 = str2+":"+str2;
    str3 = str3+":"+str3;
    str4 = str4+":"+str4;
    this->selec(str1, str2, M1, str3, str4);

}


// Não funciona direito porque toda vida que é finalizado o método a função está entendendo que deve apagar o parÂmetro de entrada.
template <class TipoClasse>
void matrix<TipoClasse>::operator=(matrix M1)
{
//    if(this->lin == 0 && this->col == 0)
    this->init(M1.lin, M1.col);
    try{
        for(  int i = 0; i < this->lin; i++)
            for(  int j = 0; j < this->col; j++)
                 this->Mat[i][j] = M1.Mat[i][j];
        }
    catch(int i)
    {
        cout << "não é possível realizar a atribuição"<< endl;
    }
}

template <class TipoClasse>
void matrix<TipoClasse>::operator= (string valor)
{
    this->init(valor);
}

//-----------------------------------------------------------------//
//-----------------------------------------------------------------//
//-------------------Mostra as matrizes em tela -------------------//
//-----------------------------------------------------------------//
//-----------------------------------------------------------------//

template <class TipoClasse>
void matrix<TipoClasse>::print()
{
    cout << endl;
    for(  int i = 0; i < this->lin; i++)
    {
        for(  int j = 0; j < this->col; j++)

             cout << this->Mat[i][j] <<"  ";

         cout <<'\n';
    }
}

#endif // MATRIX_H
